Simulations of Two-Dimensional Melting on the Surface of a Sphere 
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We have simulated a system of classical particles confined on the surface of a sphere interacting with 
a repulsive r -12 potential. The same system simulated on a plane with periodic boundary conditions 
has van der Waals loops in pressure-density plots which are usually interpreted as evidence for a 
first order melting transition, but on the sphere such loops are absent. We also investigated the 
structure factor and from the width of the first peak as a function of density we can show that the 
growth of the correlation length is consistent with KTHNY theory. This suggests that simulations 
of two dimensional melting phenomena are best performed on the surface of a sphere. 
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Halperin, Nelson jjj and Young |2) established several 
years ago a theory of defect-mediated melting of two di- 
mensional (2D) crystals. It is based on ideas of Kosterlitz 
and Thouless || . It is envisaged that the transition from 
crystal to liquid takes place via two continuous transi- 
tions instead of the single first order melting transition 
of three dimensional crystals. A hexatic phase appears 
between ordered (solid) and isotropic (liquid) phases. We 
shall refer to this theory as KTHNY theory. The crys- 
talline phase has both long-range crystalline order and 
bond orientational order, but the crystalline order is lost 
at the transition to the hexatic phase when according to 
KTHNY theory dislocation pairs unbind. As the tem- 
perature is raised further a second continuous transition 
occurs above which orientational order is lost when discli- 
nation pairs unbind and the isotropic liquid state forms. 
The correlation length £ over which there is short-range 
crystalline order diverges as T approaches the melting 
temperature. It is given by the expression: 



£(T) oc exp 



(1) 



{T/T m -l)\ 

where & is a numerical factor and v = 0.36963 . . . 

The KTHNY scenario is not the only mechanism pos- 
sible for two-dimensional melting. It may happen that 
a single first order transition occurs before the defects 
unbind. Experimentally, there are systems to which the 
KTHNY theory applies and there are systems where just 
a single first order transition takes place. Electrons on 
the surface of helium p| , submicron polymer colloids con- 
fined between glass plates || , structural order of two di- 
mensional charge-density-waves in Nb x Tai— X S2 J(| and 
colloidal particles confined to a monolayer with dipole 
interactions |7| are examples of the former while xenon 
on graphite His a example of the latter. 

However, Monte Carlo (MC) and molecular dynam- 
ics (MD) simulations of two dimensional melting almost 
always seem to indicate that melting takes place via a 
first order transition |9"Hl2|. Recently, there has been a 
very large scale simulation |i"3| of a system of particles 
interacting via a r -12 potential with periodic boundary 



conditions. Systems of 4096, 16384 and 65536 particles 
were studied and as the number of particles increased 
the size of the van der Waals loops decreased implying 
that in the thermodynamic limit the melting transition 
might be continuous. It is the chief purpose of this paper 
to point out that it is much easier to get results valid in 
the thermodynamic limit by making the two-dimensional 
system live on the surface of a sphere. 

Problems arise with MC and MD simulations when the 
relaxation times in the system become comparable with 
the timescale of the simulation. When plotting pressure 
versus density along an isotherm, a hysteresis loop can 
appear around the phase transition region. This is ex- 
pected for systems with a first order transition but a sys- 
tem out of equilibrium can also show hysteresis. It has 
been found that deliberately introducing vacancies into 
the simulations with periodic boundary conditions de- 
creases the jump at the first order transition E3. This 
suggests that one of the problems of doing simulations 
with periodic boundary conditions is that the timescales 
employed in these simulations might be less than the 
timescale for nucleating defects such as vacancies, dis- 
location pairs etc. thus preventing the attainment of 
true equilibrium. Now on the surface of a sphere the 
crystalline state always contains at least 12 disclinations 
(i.e 5-fold coordinated sites) due to Euler's theorem (see 
eg. Cq|). Furthermore, dislocations appear even in the 
ground state to screen the 12 disclinations in order to re- 
duce their elastic strain energy fl5fl . In other words, the 
topology of the sphere always forces a number of defects 
into the crystalline state. But instead of being a dis- 
advantage as one might have first thought, the presence 
of the defects seems to allow the system to relax more 
readily and so explore more phase space on the timescale 
of the simulation. It is possible that this is the mecha- 
nism which makes simulations on the surface of a sphere 
closer to the thermodynamic limit than those done with 
periodic boundary conditions 

Thermodynamic properties such as the free energy per 
particle for particles with short range interactions are 
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the same on the flat plane and the sphere as N — > oo. 
When N — > oo the local curvature vanishes (R, the 
radius of the sphere is varied with TV so as to keep the 
surface density p constant) so that locally the system on 
the sphere appears flat. 

We have studied by means of MD a system of N par- 
ticles constrained to move on the surface of a sphere. 
Particles interact with each other by the repulsive pair 
potential 



v(rij) = e{a/r i: j) 12 



(2) 



The distance between particles i,j is measured along 
the geodesic truncating at r%j = 3a. m is given by the 
expression 



T{j 17 ij R 



(3) 



where i9y is the angle between particles i and j. 

■dij — arccos (sin 6i sin 6j cos (fy — <j)j) + cos 0+ cos 6j ) (4) 

In this paper we will use reduced units (a = e = m = 1). 

We employed a velocity Verlet algorithm. It was 
adapted to our particular case of polar coordinates 9i, & 
with < 0i < 7T and < 6i < 2n for i = 1, 2, . . . N: 



i(t + St) = 0i(t) + {v 6i (t)St + l/2a gi (t)St 2 )R- 



(5) 



YK ' Ynj sin0i(t + St) + sin 0i(t) 

where a ( f >i — v ( f H ,ag i = vg i ,v ( p i = Rsm.0i4>i,v$ i = R6i 
Accelerations , ag i for this problem are as follows 
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Temperature is introduced in the simulation by rese- 
lecting the velocities of all the particles at once according 
to a Boltzmann distribution. This reselection is done at 
equally spaced intervals of time . 

It seems natural to choose values for N such that the 
particles can be disposed on the sphere in a triangular- 
like lattice. This occurs when N = 10T + 2 where T is 
equal to 



T = h 2 + hk + k 2 



(7) 



with h and k integers |17| , |18| . For values of TV's satisfing 
Eq. (^) the particles can be arranged on the surface of 
the sphere into a state with full icosahedral symmetry. In 
this paper we study systems of 72, 122 and 272 particles. 



It easy to see that the calculation of the accelerations 
from Eq. (^) is slower than the corresponding calcula- 
tions on a flat plane as each element in the sum in (||) 
requires the (slow) computation of several trigonometric 
functions. This limits our simulation to rather small val- 
ues of N, so it is fortunate that on the sphere results 
consistent with KTHNY theory can be seen with small 
values of N. The obvious idea of using look-up tables 
for the trigonometric functions led to considerable errors. 
We expect though that improvements in our numerical 
procedures can be found. This would allow us to study 
larger systems and so see the hexatic phase, which be- 
cause it exists only in a narrow density region above the 
melting density, is hard to disentangle from finite size 
effects in small systems. 
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FIG. 1. Pressure-density plot for 72 (solid circles), 122 
(empty circles) and 272 (empty squares) particles on a sphere. 
Results from Broughton et al. (crosses) for the same potential 
on a flat plane are also plotted for comparison. 

We have calculated the pressure P and the structure 
factor along the isotherm T = 1. A simulation time 
400,000(W was used for each value of the density, where 
St = 0.005(mcr 2 /e) 1 / 2 . The pressure was evaluated using 
the expression |19[ 



G 



N 



i<j 



(8) 



Fig. 1 shows the pressure-density isotherm for 72 (solid 
circles), 122 (empty circles) and 272 (empty squares). 
Crosses represent data obtained by Broughton et al. M 
for the same interacting potential on a flat plane. At low 
densities, the pressures are the same for the flat plane and 
for the sphere when both systems are closer to ideal gas 
behavior. As one can see in Fig. 1, there is no hysteresis 
loop on the sphere indicating that the solid phase does 
not melt by a first order transition. 

If one concedes that there is no evidence for a first order 
melting transition then one might try to argue that what 
is happenning instead is that the topology of the sphere 
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is preventing a solid phase forming and that the system is 
always liquid. However, by studying the structure factor 
we can see that a transition to a crystalline phase does 
indeed occur on the sphere. 

The structure factor S(k) is the Fourier transform of 
the pair correlation function h(r). h(r) is related to the 
pair distribution function g(r) simply by h(r) = g(r) — 
1. Several authors have developed methods to calculate 
the structure factor in spherical geometries p(i| , pT| |. We 
obtained it from the pair distribution function using the 
equation 

/•7T 

S(k) = 1 + 2npR 2 / h{R0) sm9J {kR9)dd (9) 
Jo 

Jo is the Bessel function of zeroth order. Fig. 2 shows 
the structure factor for p — .9. The dashed line is a 
Lorentzian fit to the first peak of the structure factor 
(which occurs at a wave- vector corresponding to the first 
reciprocal lattice vector of a triangular lattice). We define 
the inverse correlation length as the width of the first 
peak. Fig. 3 shows the behavior of the inverse width as 
a function of the density. 



0.8 r— 

0.7 - 

0.6 - 

0.5 - 

0.4 - 

0.3 - 

0.2 - 

0.1 - 

0.0 — 
0.8 



*!-. 



0.9 1.0 



1.1 1.2 



1.3 



FIG. 3. Width of the first peak of the structure factor, that 
is the inverse correlation length, as a function of density p. 
N is equal to 72 (empty circles), 122 (solid squares) and 272 
(solid circles) 

£ saturates when p w 1 when £ is of the order of the ra- 
dius R of the sphere. This kind of behavior is as expected 
for a continuous phase transition. 

Eq. (|l|) gives the dependence of £ in the liquid phase 
along an isochore. The dependence of £ along an isotherm 
is given by 
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FIG. 2. Structure factor for p = 0.9 (solid circles) . Dashed 
line is the Lorentzian fit used to obtain the width of the 
first peak, a is a measure of the particle separation where 
a = (np)' 



-1/2 



£(p) oc exp 



{{Pm/pf - 1) 



(10) 



where v = 0.36963 as before and p m is the melting den- 
sity. 

To obtain Eq. (|o|) from Eq. ([l]) we have used a 
scaling properties of the r~ 12 potential Q; viz that sys- 
tems with different T's and p's share the same thermody- 
namic properties if they have the same value of T, where 
r = (irp) 6 /(k B T). 

Fig. 4 shows a log-log plot of ln£ as a function of 
(1/p 6 - l) - 36963 . (We have taken p m = 1). A straight 
line of slope equal to -1 is also shown. This slope is the 
KTHNY prediction. The plot should thus should have 
a slope of -1 for densities approaching p m from below. 
(Very close to p m finite size effects will modify the be- 
havior). The results obtained are clearly consistent with 
KTHNY theory, Eq. @. 
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FIG. 4. ln£ versus (l/p B - 1)0-36963 The glope of the 
straight line is -1 which corresponds to KTHNY predictions. 
N is equal to 72 (empty circles), 122 (solid squares) and 272 
(solid circles) 

To summarize, we have simulated a system of classical 
interacting particles on a sphere interacting via a short- 
range interaction. In this geometry there is no hysteresis 
loop as on a flat plane. By studying the structure factor, 
evidence of a continuous melting transition was found 
near p = 1. The correlation length of short-range crys- 
talline order in the liquid phase diverged on approach- 
ing the transition as predicted by KTHNY theory. We 
believe that this constitutes strong evidence that simu- 
lations of two dimensional melting phenomena are best 
performed on the surface of a sphere and implies that the 
first order transition so often reported in simulations on 
the flat plane is nothing but a numerical artifact. 
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